#include <iostream>
#include <vector>
#include <string>
#include <math.h>
#include <lapacke.h>

using namespace std;

struct Point // 定义点的结构，存储坐标信息
{
    double x;
    double y;
    Point():x(0),y(0){};
    Point(double x, double y):x(x),y(y){};
};

struct Value // 存储函数信息
{
    Value(){};
    double F0(Point g); // u的原函数
    double F1(Point g); // u的一阶导函数
    double F2(Point g); // u的二阶导函数
};

class Cricle // 圆
{
public:
    Point O; // 圆心
    double R; // 半径
    Cricle():R(0),O(0,0){};
    Cricle(Point _O, double _R)
    {
	O = _O;
	R = _R;
    };
    bool InOnCricle(Point P);// 判断点P是否在圆内或圆上
    vector<double> Calculate_X(double Y); // 给定Y坐标求圆对应的两个X坐标
    vector<double> Calculate_Y(double X); // 给定X坐标求圆对应的两个Y坐标
};

class Grid
{
private:
    int N;
    Cricle C;
public:
    Grid()
    {
	N = 0;
	Cricle _C;
	C = _C;
    };
    Grid(int _N, Cricle _C)
    {
	N = _N;
	C = _C;
    };
    vector<Point> G,IG,IJ;
    vector<vector<Point>> G0;
    vector<vector<double>> Label;
    Point get(int i, int j);
    void get_RDG(); // 获得规则的所有点的信息
    void get_IDG(); // 获得需要计算点的信息，并给所有点打上标记
};

class Solver
{
private:
    int N;
    string TYPE;
    double R,x,y;
public:
    Solver()
    {
	R = 0;
	x = 0;
	y = 0;
	N = 0;
	TYPE = "D";
    }
    Solver(int _N,string _TYPE)
    {
	R = 0;
	x = 0;
	y = 0;
	N = _N;
	TYPE = _TYPE;
    };
    Solver(int _N,string _TYPE, double _x, double _y, double r)
    {
	R = r;
	x = _x;
	y = _y;
	N = _N;
	TYPE = _TYPE;
    };
    Point O;
    Cricle C;
    Grid grid;
    vector<double> e;
    vector<double> F,U,F0;
    vector<vector<double>> A;
    void get_AF();
    void get_U();
    void get_F0();
    void ERROR();
    void ERROR1();
    void solver()
	{
	    if (TYPE == "D")
	    {}
	    else if (TYPE == "N")
	    {
		cout << "Sorry. Can't caculate Neumann." << endl;
		exit(0);
	    }
	    else if (TYPE == "M")
	    {
		cout << "Sorry. Can't caculate Mixed." << endl;
		exit(0);
	    }
	    else
	    {
		cout << "Error boundary condition type input." << endl;
		exit(0);
	    }
	    Point _O(x,y);
	    O = _O;
	    Cricle _C(O,R);
	    C = _C;
	    Grid _grid(N,C);
	    grid = _grid;
	    grid.get_IDG();
	    get_AF();
	    get_U();
	    get_F0();
	    ERROR();
	};
};
//---Value----------------------------------------------
double Value :: F0(Point g)
{
    double f;
    f = exp(g.y + sin(g.x));
    return f;
};

double Value :: F1(Point g)
{
    double f;
    f = exp(g.y + sin(g.x)) + cos(g.x) * exp(g.y + sin(g.x));
    return f;
};

double Value :: F2(Point g)
{
    double f;
    f = exp(g.y + sin(g.x)) * (sin(g.x) * sin(g.x) + sin(g.x) - 2);
    return f;
};
//---Cricle---------------------------------------------
bool Cricle :: InOnCricle(Point P)
{
    if ((O.x-P.x)*(O.x-P.x)+(O.y-P.y)*(O.y-P.y) <= R*R)
	return true;
    return false;
};

vector<double> Cricle :: Calculate_X(double Y)
{
    vector<double> X;
    X.resize(2);
    X[0] = O.x - sqrt(R*R - (O.y - Y)*(O.y - Y));
    X[1] = O.x + sqrt(R*R - (O.y - Y)*(O.y - Y));
    return X;
};

vector<double> Cricle :: Calculate_Y(double X)
{
    vector<double> Y;
    Y.resize(2);
    Y[0] = O.y - sqrt(R*R - (O.x - X)*(O.x - X));
    Y[1] = O.y + sqrt(R*R - (O.x - X)*(O.x - X));
    return Y;
};
//---Grid-----------------------------------------------
Point Grid :: get(int i, int j)
{
    Point g;
    g.x = (double)i/N;
    g.y = (double)j/N;
    return g;
};

void Grid :: get_RDG()
{
    G.resize((N + 1)*(N + 1));
    G0.resize(N + 1);
    for (int i = 0; i < G.size(); i++)
    {
	G[i] = get(i/(N + 1), i%(N + 1));
    }
    for (int i = 0; i < G0.size(); i++)
    {
	G0[i].resize(N + 1);
	for (int j = 0; j < G0[i].size(); j++)
	{
	    G0[i][j] = get(i,j);
	}
    }
};

void Grid :: get_IDG()
{
    get_RDG();
    Label.resize(N + 1);
    Point ij;
    int j = 1;
    for (int i = 0; i < Label.size(); i++)
    {
	Label[i].resize(N + 1);
    }
    for (int i = 0; i < Label.size()*Label.size(); i++)
    {
	if (G[i].x != 0 && G[i].y != 0)
	    {
		if (G[i].x != 1 && G[i].y != 1)
		{
		    if (!C.InOnCricle(G[i])) // 如果点在圆内标记为-1
		    {
			ij.x = (N)*G[i].x;
			ij.y = (N)*G[i].y;
			IJ.push_back(ij); // 存储第j个需要计算的点坐标
			IG.push_back(G[i]);// 存储第j个需要计算的点
			Label[ij.x][ij.y] = j;
			j++;
		    }
		    else // 点需要计算，标记为j
		    {
		        Label[(N+1)*G[i].x][(N+1)*G[i].y] = -1;
		    }
		}
	    }
    }
    if (j != (N-1)*(N-1)+1 && j > (N-1)*(N-1)-3)
    {
	cout << "Error cricle input." << endl;
	exit(0);
    }
};
//---Solver---------------------------------------------
void Solver :: get_F0()
{
    Value V;
    F0.resize(grid.IG.size());
    for (int i = 0; i < F0.size(); i++)
    {
	F0[i] = V.F0(grid.IG[i]);
    }
};
void Solver :: get_AF()
{
    vector<Point> IG,IJ;
    vector<vector<Point>> G0;
    vector<vector<double>> Label;
    IG = grid.IG;
    IJ = grid.IJ;
    Label = grid.Label;
    G0 = grid.G0;
    Value V;
    int m = IG.size();
    F.resize(m);
    A.resize(m);
    for (int i = 0; i < m; i++)
    {
	F[i] = -V.F2(IG[i]);
	A[i].resize(m);
    }
    for (int i = 0; i < m; i++)
    {
	vector<Point> theta;
	double u = 0;
	theta.resize(2);
	for (int k = 0; k < 2; k++)
	{
	    theta[k].x = 0;
	    theta[k].y = 1;
	}
	// 探测上方的点
        if (Label[IJ[i].x - 1][IJ[i].y] == 0)
	{
	    u = V.F0(G0[IJ[i].x - 1][IJ[i].y]);
	}
	else if (Label[IJ[i].x - 1][IJ[i].y] == -1)
	{
	    Point u;
	    theta[0].y =N*(IG[i].x - C.Calculate_X(IG[i].y)[1]);
	    u.x = C.Calculate_X(IG[i].y)[1];
	    u.y = IG[i].y;
	    theta[0].x = V.F0(u);
	    F[i] = F[i] + 2*N*N*theta[0].x/(theta[0].y*(1+theta[0].y));
	}
	else
	{
	    A[i][Label[IJ[i].x - 1][IJ[i].y] - 1] = -2*N*N/(1+theta[0].y);
	}//
	// 探测下方的点
	if (Label[IJ[i].x + 1][IJ[i].y] == 0)
	{
	    u = V.F0(G0[IJ[i].x + 1][IJ[i].y]);
	}
	else if (Label[IJ[i].x + 1][IJ[i].y] == -1)
	{
	    Point u;
	    theta[0].y =N*(- IG[i].x + C.Calculate_X(IG[i].y)[0]);
	    u.x = C.Calculate_X(IG[i].y)[0];
	    u.y = IG[i].y;
	    theta[0].x = V.F0(u);
	    F[i] = F[i] + 2*N*N*theta[0].x/(theta[0].y*(1+theta[0].y));
	}
	else
	{
	    A[i][Label[IJ[i].x + 1][IJ[i].y] - 1] = -2*N*N/(1+theta[0].y);
	}//
	F[i] = F[i] + 2*N*N*u/(1+theta[0].y);
	u = 0;
	// 探测左侧的点
	if (Label[IJ[i].x][IJ[i].y - 1] == 0)
	{
	    u = V.F0(G0[IJ[i].x][IJ[i].y - 1]);
	}
	else if (Label[IJ[i].x][IJ[i].y - 1] == -1)
	{
	    Point u;
	    theta[1].y =N*(IG[i].y - C.Calculate_Y(IG[i].x)[1]);
	    u.y = C.Calculate_Y(IG[i].x)[1];
	    u.x = IG[i].x;
	    theta[1].x = V.F0(u);
	    F[i] = F[i] + 2*N*N*theta[1].x/(theta[1].y*(1+theta[1].y));
	}
	else
	{
	    A[i][Label[IJ[i].x][IJ[i].y - 1] - 1] = -2*N*N/(1+theta[1].y);
	}//
	// 探测右侧的点
	if (Label[IJ[i].x][IJ[i].y + 1] == 0)//右
	{
	    u = V.F0(G0[IJ[i].x][IJ[i].y + 1]);
	}
	else if (Label[IJ[i].x][IJ[i].y + 1] == -1)
	{
	    Point u;
	    theta[1].y =N*(- IG[i].y + C.Calculate_Y(IG[i].x)[0]);
	    u.y = C.Calculate_Y(IG[i].x)[0];
	    u.x = IG[i].x;
	    theta[1].x = V.F0(u);
	    F[i] = F[i] + 2*N*N*theta[1].x/(theta[1].y*(1+theta[1].y));
	}
	else
	{
	    A[i][Label[IJ[i].x][IJ[i].y + 1] - 1] = -2*N*N/(1+theta[1].y);
	}//
	F[i] = F[i] + 2*N*N*u/(1+theta[1].y);
	A[i][i] = 2*N*N*(1/theta[0].y + 1/theta[1].y);
    }
};

void Solver :: get_U()
{
    U = F;
    vector<double> A0;
    for (int i = 0; i < A.size(); i++)
	{
	    for (int j = 0; j < A[i].size(); j++)
	    {
		A0.push_back(A[j][i]);
	    }
	}
	lapack_int info, m, n, lda, ldb, nrhs;
	m = U.size();
	n = U.size();
	nrhs = 1;
	lda = U.size();
	ldb = U.size();

	info = LAPACKE_dgels(LAPACK_COL_MAJOR, 'N', m, n, nrhs, & *A0.begin(), lda, & *U.begin(), ldb);
};

void Solver :: ERROR()
{
    e.resize(3);
    double t;
    for (int i = 0; i < F0.size(); i++)
    {
	t = F0[i] - U[i];
	if (t < 0)
	{
	    t = -t;
	}
	e[0] = e[0] + t/(N*N);
	if (t > e[2])
	    e[2] = t;
	e[1] = e[1] + t*t/(N*N);
    }
    e[1] = sqrt(e[1]);
    cout << "n = " << N << "," << endl;
    cout << "norm_1 of error is " << e[0] << endl;
    cout << "norm_2 of error is " << e[1] << endl;
    cout << "norm_infty of error is " << e[2] << endl;
    cout << endl;
};

void Solver :: ERROR1()
{
    e.resize(3);
    double t;
    for (int i = 0; i < F0.size(); i++)
    {
	t = F0[i] - U[i];
	if (t < 0)
	{
	    t = -t;
	}
	e[0] = e[0] + t/(N*N*N*N);
	if (t > e[2])
	    e[2] = t;
	e[1] = e[1] + t*t/(N*N*N*N);
    }
    e[1] = sqrt(e[1]);
    cout << "n = " << N << "," << endl;
    cout << "norm_1 of error is " << e[0] << endl;
    cout << "norm_2 of error is " << e[1] << endl;
    cout << "norm_infty of error is " << e[2]/(N*N) << endl;
    cout << endl;
};
//------------------------------------------------------
